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Summary 

Some path-following techniques are described and 
compared with other methods. Use of multipurpose 
techniques that can be used at more than one stage 
of the path-following computation results in a system 
that is relatively simple to understand, program, and 
use. Comparison of path-following methods with the 
method of parametric differentiation reveals definite 
advantages for the path-following methods. The fact 
that parametric differentiation has found a broader 
range of applications indicates that path-following 
methods have been underutilized. 

Introduction 

Considerable literature exists on path-following 
theory and methods (see, for example, ref. 1 and refs, 
included therein). An outstanding example is the 
paper by Keller (ref. 1, pp. 359-384) which includes, 
not a single metric in the domain space, but a 
general metric form, and no less than four techniques 
for locating bifurcating branches. Furthermore, the 
arc length normalization technique used by Keller 
eliminates the singularity at a turning point. 

The path-following methods used in the present 
investigation are somewhat different from those de- 
scribed by Keller. The approach is simple in ap- 
plication, yet versatile, and requires a conservative 
amount of programming. 

The great variety of problems to which some path- 
following techniques are applied preclude a categor- 
ical assertion that any one technique or system of 
techniques is optimal. A certain class of problems, 
for example, is concerned with bifurcation phenom- 
ena as such and the problems associated with obtain- 
ing convergent solutions in the vicinity of the singular 
point. For such problems, few path- following con- 
cepts are involved, and the methods described herein 
may not be helpful. These methods are directed to- 
ward problems, such as the majority of those dis- 
cussed in the references, that involve path-following 
techniques: predictor and corrector methods; loca- 
tion and classification of singular points; methods 
for stepping through singular points; and locating 
secondary branches and disconnected solution sets. 

The path-following literature appears to have fo- 
cussed on the bifurcation theory and a restricted 
class of physical problems that have bifurcating so- 
lutions, such as beam buckling (Reiss, ref. 1, pp. 37- 
71), chemical reactors (Ray, ref. 1, pp. 285-315), 
Benard convection (ref. 2, pp. 562-565), and Taylor- 
Couette flow (ref. 3, pp. 575-581). On the other 
hand, the closely related method of parametric dif- 
ferentiation (MPD) has been applied to a greater 
variety of problems. A partial listing of a 1985 li- 


brary computer search revealed over 30 publications 
reporting on MPD analysis of problems involving ra- 
diation gas dynamics, ignition and combustion, so- 
lution of nonlinear Volterra integro-differential equa- 
tions, nonlinear noise propagation, transonic airfoil 
and wing flows, unsteady transonic flows, internal 
transonic flow, compressible boundary layers, bound- 
ary layers with blowing, magnetofluid dynamics, and 
structural optimization. 

However, there are inherent problems with MPD, 
both theoretical (it yields an approximate solution) 
and in its implementation. It will be shown that 
these problems are overcome by using path-following 
methods. Consequently, a large class of problems 
becomes subject to essentially exact solution. 
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real Hilbert space 

value of determinant obtained by 
discretizing L 

function defining nonlinear source 
term in equation (16a) 

function of x 

— 6f, variation of / 

linear operator obtained by 
linearization of Q 

method of parametric 
differentiation 


P, G, z , q , H functions used to define F 

(eqs. (16b) through (16g)) 

Q differentiable nonlinear operator 

onBxi? 


P, R n 
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A 

6 
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Ac 


real one- and n-dimensional space 

perturbation function 

“arc length” on solution path in 
B x R 

independent variable 
increment in parameter 
variation of function or operator 
coefficient (eq. (24)) 
homotopy parameter 
real parameter 

values of A that bracket A c , A a < 
< A b 

estimated or actual critical 
(singular point) value of A 


residual of nonlinear operator 
equation 

denotes the L2" norm m real 
Hilbert space 

number of step along path 

starting condition 

conditions after first and second 
calculation steps, except in 
equations (29) through (31) 
where the subscripts denote 
different parameters 

Primes with a symbol indicate differentiation with 
respect to independent variable x. 

Analysis 

Basic Path-Following Equations 

It is assumed that the nonlinear operator Q[f, A] 
is defined on B X R where B is a real Hilbert space 
and R is the real line and that Q is continuously 
differentiable with respect to A and has a Gateaux 
variation with respect to /. It is required to solve 


P 

ii ii 

Subscripts: 

i 

0 

1,2 


Since the first term on the right vanishes, equa- 
tion (1) is satisfied approximately for /i, Ai provided 
that 

£[/o,Mo] + ^AA = 0 ( 8 ) 

Equation (8) is a linear equation for the variation 
h which, by equation (5), yields the approximate 
solution fi at Ai = Aq + A A. If this process is 
continued several times until A is no longer close to 
Aq, then the cumulative error becomes unacceptable. 
Therefore, it is necessary to provide a means of 
eliminating the error at any step along the solution 
path. This elimination is accomplished by a corrector 
calculation as follows. 

Calculate the residual of the nonlinear equa- 
tion (1): 

p(/l,^l) = Q [/li Ai] ( 9 ) 

Hold A constant at A = Ai and calculate the variation 
h that satisfies 


L[fi,h;\i] - -p 

Then the function f\ c = f\ + h satisfies 


(10) 


Q [fie, Ai] « Q [/I, Ai] + L [/, , h\ Ai] - 0 (11) 


Q[f, A] =0 (1) 

for a range of values of A, given a starting solution 
fo{x) at Aq: 

Q[/o, A 0 ] = 0 (2) 

(It will be shown later that it is sufficient to have only 
an approximate solution to begin the procedure.) If 
the parameter is incremented to Ai = Ao + AA then, 
in general, 

Q[/o,Ai]«Q[/o,A 0 ] + ^AA^O (3) 

The first variation of the operator Q\ 

6Q = L[f 0 ,h ;\ 0 ] (4) 

is a linear operator operating on the variation h(x) 

of f{x), 

fl = fo + h (5) 

Then 

Q [/i> A 0 ] « Q [/o, A 0 ] + L [/o, h; A 0 ] (6) 

If both / and A are incremented simultaneously, then 


This procedure can be iterated to drive the residual 
to machine zero at A = X\. 

It is important to observe that neither equa- 
tion (8) nor equation (10) need be formulated in 
terms of the inverse operator L -1 . They can be 
solved by any convenient numerical formulation, ei- 
ther explicit or implicit. Thus, near-singular prob- 
lems, in which the highest derivative term is small, 
can be treated without the difficulties that arise in 
such methods as that of Adomian (ref. 3). 

Singular points occur where the linear operator L 
vanishes: 

L[f,h; A] = 0 (12) 

A singular point may correspond to a bifurcation 
or a turning point of the solution path in the 
cross-product domain space. Considerable literature 
(ref. 4, for example) exists on the mathematics of bi- 
furcation and solutions in the vicinity of a singular 
point. This theory is not repeated herein inasmuch 
as the primary concern of the present study is that 
of specific path-following computational techniques. 

It is useful to define a kind of arc length to 
describe distances along the solution path in B X R. 
Let || || denote the L 2 -norm in B , 


Q [/i. All » Q [/o, Ao] + L [fo, h- Ao] + H AA ( 7 ) 



( 13 ) 
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Then the arc length in the domain space B x R can 
be defined by 

ds 2 = P/ll 2 + (dX) 2 (14) 

If several parameters Ai, A 2 , A n are involved, 
the ordered array of parameters represents a vector 
A in R n . Then the domain space is B x R n and the 
arc length is 

ds 2 = \\6ff + \d\\ 2 (15) 

Equation (14) and its approximation obtained by re- 
placing the differentials with finite increments are 
somewhat similar to the arc length normalization 
equations (2.10c) and (2.10d) of Keller (ref. 1, 
pp. 359-384). Setting 0 = ^ in equation (2.10c) of 
Keller yields one half the value of ds 2 given by equa- 
tion (14). 

Path-Following Procedures 

One of the primary goals of the present investiga- 
tion was to establish a set of procedures that would 
be (1) relatively simple to understand, (2) simple 
to describe mathematically, and (3) require a fru- 
gal amount of computer code. Although it is recog- 
nized that some relatively recalcitrant problems may 
require more complex procedures, the methods de- 
scribed here have proved effective for the problems 
studied. 

Model problem . To facilitate understanding of 
the methods to be described, the procedures will be 
illustrated by referring to a specific model problem. 
This example is of the same general type as that 
described by Keller (ref. 1, pp. 359-384), but with 


a different nonlinear term 
d 2 f 

-^ + F(f,X,x) = 0 (16a) 

F = 2q(X) + n 2 H(X) P{G[z(f,X,x)}} (16b) 

P(G) = sin G (16c) 

G(z) — z - z 2 (16d) 

z{f,X,x) = f - q{X) x( 1-x) (16e) 

q (X) = A 2 e“V2 (I6f) 

H{ X) = X (16g) 

/(0) = /(l) = 0 (I6h) 


For this problem, the predictor equation (8) is 

<Ph d£dPdG (Wdq dF dH\ A X 

dx 2 + dPdG dz h + V dq dX + dH dX ) ^ ~ 0 

(17) 

In equation (17), the first two terms represent the 
linear first variation of the operator of equation (16a) 

M/o,Mo] = !^ + §^§(*(/o.Ao))ft (18) 

The residual, after the first predictor calculation, is 
obtained from the nonlinear equation (16a): 

p(h,Xi,x) = ^ + F(F l ,X u x) (19) 

Then, the linear corrector equation is obtained by 
combining equations (18) and (19) 

L[fi,h;\\] + p(/i, Ai,x) = 0 (20) 

which is the basic equation for the Newton procedure. 

Path following in nonsingular regions . Special 
means are required to compute the solution paths 
in the vicinity of singular points, but for basic path 
following in nonsingular regions (see fig. 1), the 
predictor-corrector procedure described previously is 
effective. To complete one step along the path, the 
predictor equation (17) (fig. 2(a)) must be solved 
once and the corrector equation (20) perhaps several 
times. 

Simpler versions of this procedure have also been 
used effectively. One of these is to “omit” the 
predictor step. Beginning with the initial solution 
/o, Ao, one can increment A and compute the residual 
p(/ 0 j Ai,ar). Then, starting with arguments /o,Ai, 
iterate the corrector equation 

L[f,h;X l }+ p(f,X lt x) = 0 (21) 

until the residual is reduced to the acceptable level. 
This procedure, which appears to omit the predictor 
step, is actually equivalent to using a “horizontal” 
predictor step (see fig. 2(b)). It has been used in the 
solution of systems of algebraic equations (ref. 5). 

Another method, which provides a better predic- 
tor estimate than the horizontal predictor, but which 
still does not require the solution of a differential 
equation, is to use an extrapolation predictor. (See 
fig. 2(c).) After the first step has been completed, 
two solutions are available: /o,Aq and /i,Ai. The 
difference function fi — fo is obtained and the arc 
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length dsi is calculated by equation (14). Then A is 
incremented to A 2 , and the extrapolated predictor is 


fe — fl + 


A 2 A 

AiA 


6f 


(22) 


where 

Sf = fi - /o 

AiA = Ai - Ao 
A 2 A = A 2 — A] 

The initial solution for /o,Aq is now no longer 
required in storage and can be written over. After the 
predicted solution has been corrected to obtain the 
solution / 2 , A 2 , the second arc length element ds 2 is 
calculated. Then the increment A3A is chosen in such 
a way as to oppose any tendency for ds to increase 
or decrease, as for example 


A 3 A = A 2 A^± (23) 

This extrapolation procedure does not require ad- 
ditional programming, since the extrapolation rou- 
tine is already needed as a means for “stepping 
across” singular points, as will be explained later. 

Arc length control is incorporated automatically 
in the procedure of Keller (ref. 1, pp. 359-384), where 
the arc length normalization equation is included, 
along with the basic operator equation, as part of 
the nonlinear system to be solved. 

There exist, of course, alternate methods of per- 
forming the corrector step besides the Newton pro- 
cedure. Basic contraction mappings are discussed in 
reference 4 for solutions near a bifurcation point, and 
other procedures appropriate to specific problems are 
available. Where individual solutions of the correc- 
tor equation are costly in terms of computer time 
and storage, the corrector step is to be applied judi- 
ciously. For example, if an exact solution is required 
only at one specific value of A, several predictor steps 
may be taken and then corrected without driving the 
residual to machine zero. That is, the corrector may 
be used simply to keep the solution path from wan- 
dering too far off course until the desired value of A 
is obtained, and then the residual is driven to zero. 


convenient indicator is especially useful inasmuch as 
the vast majority of bifurcation problems of physical 
interest appear to involve simple bifurcation. 

To predict the nature of a singular point, one can 
include a subroutine that is turned on by the operator 
determinant falling below a specified level in absolute 
value. The subroutine calculates the increments 
AD in the determinant values and computes and 
stores the ratios AD /A A. If this ratio begins to 
increase in absolute value, no singularity exists. If 
it approaches zero linearly, a bifurcation is signaled. 
If it approaches zero quadratically, a turning point 
is signaled (ref. 2, pp. 490-491). Turning points 
are also signaled by the small values of AA required 
to maintain uniform arc length increments. (See 
eq. (23).) It is important to observe that these 
indicators do not determine a singularity but only 
signal a potential singular point. For example, a 
turning-point signal would occur for a solution set 
satisfying equation (12) at an inflection point. 

Computing bifurcating branches. Once a bifur- 
cation point has been detected, a number of meth- 
ods exist for locating the bifurcating branch. Keller 
(ref. 1, pp. 359-384) describes four methods. The 
first two of these methods were tried in the present 
investigation, but were eventually rejected (because 
of the complication of solving supplementary equa- 
tions) in favor of the method of perturbing the oper- 
ator to eliminate the singularity. This procedure is, 
in some respects, a simplified version of method II 
of Keller. This perturbation device is described in 
reference 2 (p. 491), and is applied by Reiss (ref. 1, 
pp. 37-71) to a buckling problem, for which the per- 
turbation has a physical significance. 

The method proceeds as follows (refer to fig. 3). 
Once a bifurcation point has been approximately lo- 
cated, say at A c , this critical value is bracketed by 
parameter values A a , A bl each several marching steps 
removed from A c , so that A a < A c < A^. Then 
the marching computation is returned to A a and a 
perturbation function is added to the operator equa- 
tion. This function satisfies the boundary conditions 
(eq. (16h)) and vanishes at A a and A5. An example 
is 


Locating and classifying singular points. Sin- 
gular points, defined as points /, A for which equa- 
tion (12) holds, are determined in the numerical al- 
gorithm by the vanishing of the determinant of the 
matrix formed by discretizing L. For simple (first- 
order) or other odd-order bifurcations, the determi- 
nant changes sign as the singularity is crossed. This 


r(e, A, x) = e (A - A a ) (\ b - A) x 2 (l - x) (24) 

Notice that the function of x chosen for the pertur- 
bation function is different from that which occurs in 
the original equation (16e). The marching is recom- 
menced at A a with a smaller step increment and with 
the perturbation function added to equation (1). If 
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the new path enters the domain of attraction of a su- 
percritical bifurcating solution, then at Aj> where the 
perturbation vanishes, the solution will be on the sec- 
ondary branch (fig. 3). If a supercritical bifurcating 
branch is not found, then the same procedure is em- 
ployed to seek a subcritical branch by starting at 
and marching backward. 

This perturbation technique possesses a theoreti- 
cal weakness. In the subspace orthogonal to the pri- 
mary path at the bifurcation point, there are an in- 
finite number of directions. The actual bifurcation 
direction may not be near the direction represented 
by the solution of the perturbed equation. When the 
perturbation function in equation (24) is prescribed 
to be in the bifurcating direction, this perturbation 
procedure is similar to Keller’s method II (ref. 1, 
pp. 359-384). 

The bifurcation direction can be obtained with 
some effort. For simple bifurcation, it is the direc- 
tion of the eigenfunction corresponding to the zero 
eigenvalue that causes L to vanish (see eq. (12)). For 
example, if the operator equation (eq. (1)) is 

f n -h A sin / = 0 
equation (12) becomes 

h n + A(cos f)h = 0 

But the bifurcation occurs from the trivial solution 
f — 0; consequently, this equation is simply 

h" + Xh = 0 

whose first nontrivial solution is at A = 7r 2 , 
h — sin nx 

Consequently, an appropriate perturbation function 
for this problem is 

r(e, 7r, z) = e (n — A a ) (A& — 7r) sin nx 

Of course, very few problems of practical interest 
can be treated analytically in this manner. How- 
ever, if a problem has been discretized so that equa- 
tion (12) appears as a matrix equation, the zero 
eigenvalue that causes the singularity can be deter- 
mined, and the bifurcation direction is that of the 
corresponding eigenvector. 

This additional work is not actually required ex- 
cept perhaps for some extremely recalcitrant prob- 
lem. It is not at all necessary for the perturbation 
function to correspond to the bifurcating direction 
but only to have some component in that direction. 


The direction is corrected by the corrector step; con- 
sequently, considerable flexibility is permissible in se- 
lecting the perturbation function. 

Another practical problem is that of selecting the 
sign and magnitude of the perturbation amplitude 
parameter e. It may be necessary to perform a trial 
calculation to determine an appropriate magnitude 
for e. Then, if the secondary branch is not found, 
a third calculation is carried out with the sign of e 
switched. 

Once a secondary branch has been located in 
this manner, one can determine if it crosses the 
primary branch or only intercepts it. If it crosses 
the primary branch, the solution can be continued 
across the primary branch by stepping across it. This 
is accomplished by prescribing a larger than usual 
extrapolation predictor when the solution approaches 
the estimated bifurcation point. Extrapolating well 
beyond the crossing point permits the corrector to 
correct to the new leg of the secondary branch (fig. 3). 

Now the bifurcation location A c can be obtained 
with good accuracy by interpolating simultaneously 
along the primary and secondary branches in the 
||/||, A plot (fig. 3). An important consideration is 
that the precise location of the bifurcation point is 
not required to locate the bifurcating branch when 
the perturbation method is used. Other methods 
(Keller, ref. 1, pp. 359-384) require the bifurcation 
point location because the computation of the bi- 
furcating branch must begin at this point. Precise 
location of the bifurcation point by a method such 
as interval-halving involves costly calculation in the 
near vicinity of the singular point where the conver- 
gence is often quite slow. 

Continuation past turning points . Two methods 
have been utilized in the present investigation to 
continue a solution path through a turning point. 
The first is to extrapolate the solution obtained 
near the estimated turning-point location. Now a 
quadratic (e.g., parabola) extrapolation is required 
rather than the linear extrapolation that is used for 
stepping through bifurcation points. (See fig. 4.) 
This procedure has been used with some success 
in the present investigation. For some problems, 
however, the extrapolations were quite large, and the 
predicted solutions were so far from the actual path 
that the corrected solution would not converge to 
the new leg of the solution path. Naturally, such an 
extrapolation will fail to yield the path continuation 
if the singular point is an inflection point. 

An alternate procedure — one that incorporates 
the perturbation technique used to locate bifurcating 
branches — has also been used successfully. Referring 
to figure 5, when the turning point is signaled, the 
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marching procedure is discontinued on the secondary 
branch and returned to the primary branch, well be- 
yond the estimated singularity location. Now, the 
marching is initiated on the primary branch toward 
the estimated singular value of A. Then, the equa- 
tion is perturbed toward the unknown leg of the sec- 
ondary branch. For the model problem (eq. (16)), 
this was accomplished by using the same perturba- 
tion function that was used to find the bifurcating 
branch but with e having a larger value and the op- 
posite sign. The perturbed solution jumps to the 
new branch (fig. 5). The perturbation is then elimi- 
nated, and the solution is continued both backward 
and forward along the new branch. 

This perturbation technique worked with a rela- 
tively small perturbation amplitude, apparently be- 
cause a secondary branch is a stronger attractor away 
from a turning point than near it. The technique is 
convenient to use since it already exists in the code 
as a means for locating a bifurcating branch. Fur- 
thermore, the perturbation method apparently is the 
only one available for locating disconnected solution 
sets (ref. 2, pp. 494-499). An interesting example of 
the occurrence of such disconnected solution sets in 
a physical context is reported in reference 4. 

The method of Keller (ref. 1, pp. 359-384)— 
eliminating the turning-point singularity by aug- 
menting the operator with the arc length normal- 
ization equation- was not used for several reasons. 
(1) It complicates the theory. (2) The singularities 
that occur for the augmented system are not nec- 
essarily the same as those for the original operator 
(ref. 1); one must monitor both systems for singular- 
ities. (3) In some problems, attempting to compute 
beyond the limiting A value leads to computational 
difficulties. For example, reference 6 reports difficul- 
ties in applying the method when the predictor step 
exceeds the limiting value of A. 

Comparison of Path Following With 
Method of Parametric Differentiation 

Path-following methods are related to several 
other procedures that attempt to solve nonlinear 
problems through a sequence of linear solutions. 
Such procedures include invariant embedding (ref. 7), 
quasi-linearization (ref. 8), the method of Adomian 
(ref. 3), and the method of parametric differentia- 
tion (MPD). Of these methods, MPD appears to have 
been used extensively. A pre-1985 library search gen- 
erated a listing of nearly three dozen reports on ap- 
plications of MPD. 

Path-following studies appear to be limited pri- 
marily to mathematical investigations of bifurcation 
and turning-point phenomena and to a restricted 
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class of physical problems that display bifurcation 
and/or turning-point behavior. The latter includes 
thermochemistry in continuously stirred chemical re- 
actors (ref. 1, pp. 285-315) and certain fluid dynam- 
ics problems of fundamental academic interest, such 
as Benard convection and Taylor-Couette flow. 

The MPD, on the other hand, has been applied 
to a large variety of problems. Some of these prob- 
lems involve singular points (e.g., ref. 9), whereas 
others use MPD to arrive at solutions that are not 
sufficiently near the starting solution for the Newton- 
Raphson procedure to converge. In some investiga- 
tions, the parameter occurs as a physical quantity 
(airfoil thickness or angle of attack, Damkohler num- 
ber, etc.), whereas in others, it is introduced arti- 
ficially. A clever use of the latter technique is de- 
scribed in reference 10. 

However, the MPD possesses several inherent 
weaknesses, the most notable of which are that it 
does not guarantee an accurate solution and it does 
not provide a means of dealing with singular points. 
One can see how these difficulties arise, and how they 
can be overcome by path-following methods, by ex- 
amining the MPD equations. Let the nonlinear op- 
erator Q depend directly and compositely on the pa- 
rameter A so that 


<3[/(A),A] = 0 (25) 


A starting solution /o,Ao is given. Differentiating 
equation (25) with respect to A yields 


L 


fo, d\ ,x 


= 0 


(26) 


where L is a linear operator operating on the function 
df/dX. The solution path is continued starting from 
/o, Aq by the relation 


= + ^ li-iAA (27) 

The similarity of equation (26) to the path-following 
equation (8) is obvious. Both are linear but the solu- 
tion of equation (8) is the variation h = 6f, whereas 
the solution of equation (26) is the derivative df/dX. 
It is seen immediately that equation (27) fails at 
singular points where df/dX is not defined. Fur- 
thermore, the error that occurs in each step is cu- 
mulative, and consequently the degree of solution 
accuracy is uncertain after the solution has been con- 
tinued for some distance. A number of attempts have 


l 

( 


i 


\ 


been made to minimize this problem. One method 
is to improve the estimate by including higher order 
terms in the series 

h = fi ~ l + d\ AX + (AA)2 + ' " (28) 

For each additional term in this series, another varia- 
tion of the operator must be taken and an additional 
(linear) differential equation must be solved (one for 
d//dA, one for d 2 //<iA 2 , etc.). 

Another method for reducing the cumulative er- 
ror is to “correct” the initial slope df fd\ at any point 
after the slope at the next point has been calcu- 
lated by some kind of averaging. Various versions 
and combinations of these methods have been used. 
Reference 11 describes four such methods, but none 
guarantees that the computed path will not wander 
away from the true solution path. Reference 11 also 
indicates that accuracy is lost if the parameter step 
length is not correlated with the numerical scheme 
step size. Reference 11 also reports some problems 
with stability of solutions. 

The path-following method attacks the problem 
of singular points by providing flags that signal their 
proximity and their nature, and by providing special 
tools for continuing the solution in the neighborhood 
of singular points. It solves the problem of solution 
accuracy by providing a corrector step for reducing 
the residual for the full nonlinear problem at any step 
along the solution path. 

Another advantage of the path-following method 
concerns problems that involve more than one pa- 
rameter. (See, for example, refs. 9 and 12.) If the 
parameters are varied simultaneously the MPD yields 


f . —L - \ . 


+ §r° u= '' 2 ’ 


•, n) (29) 


which is a set of n linear equations to be solved for 
df /d\j. The continuation equation is 


fi = h - i + Eff: aa j (30) 

The path-following method yields, for the same prob- 
lem, 


One further advantage of the path-following 
method over the MPD concerns the starting solution 
/ 0 ,A 0 . For the MPD an accurate starting solution 
is crucial. However, with the path-following method 
one can start with an approximate solution at Aq and 
apply the corrector step until an accurate solution is 
obtained at Ao- 

It is clear, therefore, that the path-following 
method has distinct advantages over the MPD, which 
has found many applications and has been applied 
innovatively in a broad range of physical and math- 
ematical problems. It appears then that the path- 
following method has been underutilized as a theo- 
retical and computational tool. 

Concluding Remarks 

Some path-following techniques have been de- 
scribed and compared with other methods. Empha- 
sis in this investigation has been on multiuse tech- 
niques that can be used at more than one stage of 
the path-following computation. The use of an ex- 
trapolation predictor, for example, both eliminates 
the requirement for solving an operator equation for 
the predictor solution and provides a means for step- 
ping across singular points. Use of the perturbation 
techniques for locating bifurcating branches reduces 
the requirement for obtaining the slowly converging 
solutions near bifurcation points. It is also useful 
for turning-point calculations and for locating dis- 
connected solutions. This incorporation of multipur- 
pose techniques results in a concise computer code 
that is relatively simple to use. 

A comparison of the method of parametric dif- 
ferentiation (MPD) with the path-following method 
indicates distinct advantages for the path-following 
method both in efficiency and in accuracy. The 
method of parametric differentiation has been ap- 
plied to a broad range of physics problems, whereas 
preoccupation with singular-point phenomena ap- 
pears to have unnecessarily restricted the path- 
following methods to the relatively narrow class of 
problems that involves singular points. 
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M/,Mi,A 2 ,--,A n ] + £|£AA,=0 (31) 

1 OA 3 

which is a single equation to be solved for the varia- 
tion h . 
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Figure 1. Structure of solution set for problem defined by equations (16a) to (16g), as indicated by plot of ||/|| 
versus A. 





(1) Solution of 
perturbed 
equation 


(2) Solution converges , 
to secondary branch ^ 
when perturbation SyZr* 
is eliminated 

s/s 

* 

Dn of /yj 

y/ 


(3) Reverse marching 
direction on 
secondary branch 


Primary 
branch - 


- Estimated 
bifurcation 
point 


(5) Precise location of bifurcation point 
by interpolating simultaneously in 
primary and secondary paths 


(4) Step across 

primary branch 
by extrapolation 


Figure 3. Perturbation approach to locating bifurcating branch and bifurcation point (after computing primary 
branch). 
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4. Method of obtaining predicted solution by polynomial extrapolation through turning point 


(3) March forward 
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Figure 5. Perturbation technique for locating remote leg of secondary branch. 
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